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Non-isotropic Attractive Bose-Einstein condensates are investigated with Newton and inverse 
Arnoldi methods. The stationary solutions of the Gross-Pitaevskii equation and their linear stability 
are computed. Bifurcation diagrams are calculated and used to find the condensate decay rates 
corresponding to macroscopic quantum tunneling, two-three body inelastic collisions and thermally 
induced collapse. 

Isotropic and non-isotropic condensates are compared. The effect of anisotropy on the bifurcation 
diagram and the decay rates is discussed. Spontaneous isotropization of the condensates is found to 
occur. The influence of isotropization on the decay rates is characterized near the critical point. 
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I. INTRODUCTION 

Experimental Bose-Einstein condensation (BEC) with 
attractive interactions was first realized in ultra cold va- 
pors of Li atoms []J , opening a new field in the study of 
macroscopic quantum phenomena. Such attractive con- 
densates are known to be metastable in spatially localized 



systems, provided that the number of condensed parti- 
cles is below a critical value Af c Q ■ More recently, Fesh- 
bach resonances in BEC of 85 Rb atoms were used to in- 
vestigate the stability and dynamics of condensates with 
two-body interactions going from repulsive to attractive 
values 0). 

Experimental atomic traps generally use a harmonic 
and slightly asymmetric potential. Thus, for most of 
the condensates produced so far, the geometry is nearly 
spherical. However, extremely asymmetric traps have 
been recently employed in experimental investigations of 
cigar-like or pancake-like condensates. 

Various physical processes compete to determine the 
lifetime of attractive condensates. The processes con- 
sidered in this paper are macroscopic quantum tunnel- 
ing (MOT) inelastic two and three body collisions 
(ICO) [lSOlIl and thermally induced collapse (TIC) 
llEl- The MQT and TIC contributions have been eval- 
uated in the literature using a variational Gaussian ap- 
proximation to the condensate wave function. However, 
this approximation is known to be in substantial quanti- 
tative error - e.g. as high as 17% for Af c H,Q,0] - when 
compared to the exact solution of the Gross-Pitaevskii 
(G-P) equation. 

In the nearly spherical isotropic case, both the elliptic 
(stable) and the hyperbolic (unstable) exact stationary 
solutions of the G-P equation were obtained numerically 
by Newton's method in ^]. These solution branches 
were shown to meet at Af c through a generic Hamiltonian 
Saddle Node (HSN) bifurcation. While the Gaussian ap- 
proximation presents an analogous HSN bifurcation, the 
amplitudes of its associated scaling laws were found to 
be in substantial (> 14%) error. A method for com- 
puting the unstable branch in the isotropic case via a 
shooting method was outlined in but generalizing 
this procedure to higher dimensions would be inefficient, 
and impossible in non-rectangular domains. The decay 
rates for the processes of MQT, ICO and TIC were also 
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computed, in the spherical case, from the numerical G- 
P solutions. They were shown to obey universal scaling 
laws. Experimentally significant quantitative differences 
were found between the exact rates and those based on 
the Gaussian approximation |l5j . 

In the extreme anisotropic cases, the variational Gaus- 
sian approximation has been computed and compared to 
the G-P solution on the elliptic (stable) branch [13, • 
This has allowed a more reliable determination of the 
critical value TVp than can be obtained by the Gaussian 
approximation 17]. However, a faithful determination 
of the lifetimes needs the computation of the hyperbolic 
(unstable) branch , which has not yet been performed 
in the anisotropic case. 

The main purpose of the present paper is to show 
that it is possible to compute the full HSN bifurcation 
diagram, and the corresponding lifetimes, in extreme 
anisotropic cases. We will do so by studying a cigar- 
like and a pancake-like condensate, and will obtain their 
MQT, ICO and TIC decay rates. While we have concen- 
trated, for simplicity, on these two axisymmetric cases, 
the new numerical methods developed in this work are 
capable of solving the general anisotropic problem. 

The paper is organized as follows. In section [H] we 
present the model considered throughout this work. Af- 
ter defining our working form of the G-P equation, we 
explain the methods that we used to obtain the station- 
ary states and their linearized stability. Section fllll is de- 
voted to the numerical determination of the bifurcation 
diagram and stability of the stationary states. Isotropic 
and non-isotropic cases are compared and the dynamics 
is discussed in terms of the HSN bifurcation. In section 
fVlwe define and compute the decay rates. Isotropic and 
non-isotropic rates are discussed, their similarity is ana- 
lyzed in terms of the spontaneous isotropization of con- 
densates. Finally, section [V] is our conclusion. Details of 
our numerical methods are given in the Appendix. 



II. PRESENTATION OF THE MODEL 
A. Gross-Pitaevskii Equation 

At low enough temperatures, neglecting the thermal 
and quantum fluctuations, a Bose condensate can be rep- 
resented by a complex wave function ^(x, t) that obeys 
the dynamics of the G-P equation 0, |2(j . Specifically, 
we consider a condensate of M particles of mass m and 
(negative) effective scattering length a in a confining har- 
monic potential V(x) = m(Cj^x 2 + Wyij 2 +Qlz 2 )/2 where 
x = (x, y, z) is the position vector. 

These variables can be rescaled with respect to any 
reference frequency Cj by using the natural quantum 
harmonic oscillator units of time to = l/u> and length 
yjh/muj. In terms of the non-dimensional vari- 
x/L , a = 4ira/L , u> x = Q x /oj, 
uj z /uj, the condensate is described 



ables t — t/r , x 
Lo y = ujy/oj and ui 2 



by the action 
A = [ dt 

with 



T = £ - fiAf 



where 



N = / d 3 x\^ 



£ = 



d 3 5 



-|W| 2 + V(x)|*| 
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The Euler-Lagrange equation corresponding to A is 
our working form of the Gross-Pitaevskii equation: 

(6) 



dt 



5* 

I V 2 - y(x) - ( |*| 



Our goal is to numerically determine the stable and un- 
stable stationary states of JJjJ and the eigenvalues of © 
linearized about these stationary states. We will carry 
out this calculation for various values of a cylindrical po- 
tential defined by UJ r — UJ X — UJy and to z : the isotropic 
case uj r = u) z , a cigar case uj r /b = lo z and a pancake case 
uj r = tjJ z jh. We will then use these results to calculate the 
condensate decay rates and compare these decay rates to 
those produced by the Gaussian approximation. 



B. Stationary States 

Stationary states of © corresponding to minima of £ 
at a given value of M can be obtained by integrating to 
relaxation the diffusion equation 



9* 
~dt 

m 
et 



8JF_ 



-V 2 -V(x)-(a|vI/| 2 - M ) 



* (7) 

(8) 



using initial data ^{t = 0) with a total number of par- 
ticles Af. The condition (JSJ fixes the value of the La- 
grange multiplier \i during the relaxation. This relax- 
ation method yields both the solution ^ and the La- 
grange multiplier [i. It is equivalent to that used in [T^] 
and can only reach the stable stationary solutions of 10. 
Unstable stationary solutions to © and Q are obtained 
by a Newton branch following method detailed in the 
Appendix. 

Note that the Lagrange multiplier /i can only affect the 
solutions of © through a homogeneous rotating phase 
factor e lfJ,t , in contrast to its particle number conservation 
effect on equation Q). However, every stationary solution 
to © is indexed by the unique /i-value that makes it 
time- independent, as shown in figure ^ 
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FIG. 1: Particle number N as a function of /i for the ex- 
act solutions (solid curves) and the Gaussian approximation 
(dashed curves) presented in section IIIII From top to bot- 
tom: pancake (u r — lj z /5), cigar (uv/5 = lu z ), and isotropic 
(uj r = Lu z ) geometries. 



C. Linearized Stability 

We now turn our attention to computing the linear 
stability of the G-P equation about a stationary state. 
We first write © in the abbreviated form: 



where 



9* 

-i— = LV + WCSf) 
at 



(9) 



(10) 



W(V) = [-V(x) - a|*| 2 + fi] * (11) 

The stationary states of © satisfy: 

= L* + VF(*) (12) 

Without loss of generality, ^ can be chosen to be real. 
Our objective is to calculate the eigenpairs (A, ip) of the 
operator that results from linearizing 10 about a sta- 
tionary state (We use \& to designate solutions to the 
nonlinear problem (|12f> and ip to designate eigenvectors, 
which are solutions to the linear problem to be defined 
below.) In order to correctly formulate the linear stability 
problem, it is necessary to first decompose ip = ip R + iip I . 
We write the linearized evolution equation 



dip 
~dt 



i(L + DW(V))(ip R : +iip I ) 



(13) 



where DW(^) is the Frechet derivative, or Jacobian, of 
W evaluated at 'J. DW(^) acts on ip via: 



DWip = DW R ip R + iDW'ip 



rl.iJ 



(14) 



where we have omitted the functional dependence of 
DW, DW R , and DW 1 on and where 



DW R = fi-V(ic)-3at> 2 
DW 1 = fi-V(x) - a* 2 . 

Equation (|13|l is then written in matrix form as 



(15a) 
(15b) 



d_ 

m 



ip R 
ip 1 



o 

DW R 



-{L + DW 1 ) 



L + DW 
The eigenmodes (\,ip R ,ip I ) satisfy 

A 



i> R 



iP R 
ip 1 



L 



-(L + DW 1 ) 
DW R 



iP R 
ip 1 



(16) 



(17) 



Note that this eigensystem is usually presented in the 
literature (se< 
the variables 



literature (see, for example, reference 22J) in terms of 



pre s 

E3) 



(co B , i/>, ip*) = (— »A, ip R + ty 1 , *p R - iip 1 ) 



(18) 



as the equivalent Bogoliubov-de Gennes coupled equa- 
tions 



UJ 



where 



<P 



L + DW B -a* 2 
a* 2 -(L-DW B ) 



<P 



/i - V(x) - 2a* 2 



(19) 



(20) 



In the following, we will work with matrix formulation 
(I17J1 because it avoids a potential notational inconsis- 
tency of 119fl arising from the fact that ip and ip* are 
complex conjugates only when u> B is imaginary. 

It is mo re c onvenient to work with the square of the 
matrix in 117(1 : 



if) 1 

(L + £W 7 )(L- 




(21) 



DW R ) 



-(L + DW R ){L + DW 1 ) 



Because (12 1 1 is block diagonal, it can be separated into 
the two problems: 



\ 2 ip R = —{L + DW r )(L + DW R )tp R 
\ 2 ip ! = —(L + DW R )(L + DW 1 )^ 1 



(22a) 
(22b) 



Problems l|22a|) and l|22b|) are closely related, Since the 
operators L, DW 1 , and DW R are all self-adjoint under 
the standard Euclidean inner product, the operators in 
(I22all and (|22b|) are adjoint to each other. If ip R is an 
eigenvector of (|22a|l with non-zero eigenvalue A 2 , then 
(L + DW R )ip R is an eigenvector of (|22b(l with the same 
eigenvalue. (Similarly, if (X,ip R ,ip I ) is an eigenmode of 
(fPT|) . then (—\,ip R ,—ip I ) is also an eigenmode of l(T7|) .) 
Thus, we solve only l(22a(l . The eigenvalues A 2 of l(22a(l - 
(|22b|l must be either complex conjugate pairs or real. We 
find them to be real and (almost all) negative, perturbed 
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only slightly from the eigenvalues of — L 2 . The eigen- 
values A of l|17l) are therefore found to be either pure 
imaginary or pure real, with most imaginary. 

Problems JTJJ and (|22a.p - l|22b[l have neutral eigen- 
modes which reflect the physical invarianccs of the prob- 
lem. Since DW 1 ^ — W(^>), then the stationary state 
$ is a neutral mode of L + DW 1 and hence of prob- 
lem J22EJ. This neutral mode is the phase mode of 10, 
since its existence is a consequence of the invariance of 
solutions \P to (|12l) under multiplication by any complex 
number on the unit circle. The corresponding eigenmode 
of problem l|22a|) is d^f/dfi. This neutral mode can be 
understood as a consequence of differentiating Jfjjl with 
respect to /1: 



0= — [(L + W)9] = (L- 



DW 1 



(23) 



Thus, 



-{L + DW r ){L + DW R )— = {L + DW 1 )-® = 0. (24) 

In terms of the original problem (|I7(I . the phase mode 
(A, ip R , ip 1 ) = (0, 0, ty) is a neutral eigenvector, while l)23|) 
shows that (A, ip R , tfj 1 ) = (0, d^f/dfi, 0) is a neutral gener- 
alized eigenvector, the two modes forming a Jordan block 
for |T7|). 

In practice, we fix \i to calculate the stationary states $ 
and the eigenvalues. The operators of (|22a[l - lj22b() depend 
on \i both explicitly and through ^ . For fi above a critical 
value fi c all eigenvalues A are imaginary, i.e. $ is an 
elliptic stationary state of ©. As n crosses fi c , we will 
see that one imaginary pair fuses at zero, and becomes 
real, with one positive and one negative value of A for 
fi < [i c . Stationary states for /1 < fi c are thus hyperbolic 
in the directions corresponding to these eigenvalues. 



III. BIFURCATION AND STABILITY OF 
CONDENSATES 
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FIG. 2: Stationary solutions of the GP equation versus 
the particle number M for the isotropic potential case with 
uo r = Q) z = a). Top: value of the energy functional £+ on the 
unstable (hyperbolic) branch and £_ on the stable (elliptic) 
branch. Bottom: square of the bifurcating eigenvalue (A±). 
Note that |A_| is the energy of small excitations around the 
stable branch. Solid lines: exact solution of the GP equation. 
Dashed lines: Gaussian approximation. 




In this section we will find the stationary solutions 
and study the stability of isotropic (uj r = u) z ), cigar- like 
(w r /5 = uj z ) and pancake-like (a;,- = lu z /5) condensates. 
These results were obtained by solving equation (|12|l for 
the stationary states and (|17|l or l)22all for the correspond- 
ing bifurcating eigenvalues. The system is discretized us- 
ing pseudo-spectral methods in a spherical domain for 
the isotropic case and in a periodic Cartesian domain for 
the non-isotropic cases. We use Newton's method to cal- 
culate the branches of stationary states. The bifurcating 
eigenvalue is found in the isotropic case by diagonalizing 
the matrix corresponding to l|17(l . In the non-isotropic 
case, we use instead the iterative inverse Arnoldi method, 
which requires only actions of the operator in (|22a|l . The 
BiCGSTAB variant of the conjugate gradient method is 
used to solve the linear systems required by both New- 
ton's method and the inverse Arnoldi method. The nu- 
merical methods we use are described in greater detail in 
the Appendix. 



FIG. 3: Condensate density \ip\ 2 as a function of radius r, 
in reduced units (see text). Solid lines: exact solution of the 
G-P equation. Dashed lines: Gaussian approximation. Stable 
(elliptic) solutions are shown for particle number J\f = 252 (a) 
and M = 1132 (b). (c) is the unstable (hyperbolic) solution 
for M = 1132 (see insert). 



A. Isotropic Condensate 

In order to compare our results with the existing ex- 
periments on quasi-isotropic condensates, we will use the 
following physical constants, corresponding to 7 Li atoms 
in a radial trap: m = 1.16 x 10~ 26 kg, a = — 27.3ao (with 
ao the Bohr radius) and Co = {u) x ujyb] z ) 1 / 3 ' = 908.41s -1 . 
These values yield a = —5.74 x 10~ 3 . With these param- 
eters, the mean-field approximation © is expected to be 
very reliable. Note that we ignore the contributions of 
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non-condensed atoms. They interact with the conden- 
sate only through a nearly constant background density 
term, inducing no significant change in the dynamics of 
the system [21j . 

The values of the energy functional £ and the (small- 
est absolute value) square eigenvalue A 2 versus particle 
number Af are shown as solid lines on figure |3 (top and 
bottom, respectively). The eigenvalues are imaginary on 
the metastable elliptic lower branch (A 2 < 0) and real on 
the unstable hyperbolic upper branch (A 2 > 0). Using 
PJl on stationary solutions we obtain d£/dAf = fi. Thus 
\x is the slope of £ and the lower branches £-, X 2 _ (re- 
spectively upper branches £+, \\) are scanned for /i > \i c 
(respectively /i < /Lt c ). The point \i = \i c determines the 
maximum number of particles Af — Af c for which station- 
ary solutions exist. We have checked that all the other 
pairs of eigenvalues are imaginary on both branches (data 
not shown). 

The dashed curves on figure [3 are derived from the 
Gaussian variational approximation which will be defined 
in section 1111 CI for the general anisotropic case. In the 
present isotropic case, this approximation can be solved 
in closed form, yielding the expressions |15| : 

4\/2^3 f-8/x + 3^/7 + 4//) 

MO*) = (25) 

7|5|(-2/i+ V7 + VJ 

£ = Af(n) (-n + 3 V7 + V) /7. (26) 



A 



B 



The number of particles Af is maximal at 7V C G = 
8v / 2^'/|5 5/4 a| for = 1/2^5. The eigenvalues 

can also be obtained in closed form from the linearized 
equations of motion |15| : 



A 2 Gu) = 8/i 2 - Any/7 + V + 2 



(27) 



By inspection of figure [21 it is apparent that both the 
solution of the GP equation and the Gaussian variational 
approximation share the same qualitative behavior, with 
quantitative discrepancies. Figure |21 shows the physical 
origin of the quantitative errors in the Gaussian approx- 
imation. It is apparent that the exact solution is well 
approximated by a Gaussian only for small Af on the 
stable (elliptic) branch. 



B. Hamiltonian Saddle Node Normal Form 

The qualitative behavior displayed on figure by the 
physical quantities £ and A 2 near the critical value Af = 
Af c is the generic signature of a HSN bifurcation defined, 
at lowest order, by the normal form p3l |24| 



q = S-(3q 2 , 



(28) 



where S = (1 —Af/Af c ) is the bifurcation parameter, j3 is 
a dimensionless constant, and q is the coordinate describ- 
ing the state of the system in the direction of the phase 
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FIG. 4: Phase portraits of the Hamiltonian saddle-node nor- 
mal form with p = q. A: S = 0.2, B: 5 = 0.1, C: S = 0. 
D: Corresponding potential $ associated to each phase por- 
trait A, B or C, with p — —d&/dq. An elliptic region bounded 
by the separatrix that starts and ends on the fixed point Q+ 
(homoclinic orbit) is present on A and B. Phase portrait C 
displays the critical merging of fixed points Q+ and Q+, and 
the disappearance of the elliptic region. 



space that becomes unstable. Indeed, introducing the 
additional non-dimensional quantities £ and 7 to define 
the appropriate energy 



1 



1 



£ = £ + -r - Sq + -/V - 7< 5, 



(29) 



2" 1 3' 

it is straightforward to derive from (|28|l that, close to the 
critical point 5 = 0, the universal scaling laws are given 
by 

£ c -£ t 6±£A6 3/2 , 



2 
± 



(30) 
(31) 



where £ c = £0, £1 = 7, £a — 2/3VJ3 and X\ = 2yf/3. 
Note that these relations can be inverted to obtain the 
parameters in (|28|l from the critical data. For the Gaus- 
sian approximation the critical amplitudes can be com- 
puted from equations l|25() and (|26|l . One finds 



£,■ 



£a = 



\2 



5 3 /> 
G4 



5 9 / 4 |a| 
4VT0. 



(32) 

(33) 
(34) 



For the exact solutions, we obtain the critical amplitudes 
£a = 1340 and X A = 14.68 by performing fits on the nu- 
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FIG. 5: Stationary solutions of the GP equation versus the 
particle number J\f for a non-isotropic potential case with 
uo r = Co and lo z — Co/5 (cigar). Top: value of the energy func- 
tional. Bottom: square of the bifurcating eigenvalue (A±). 
Solid lines: exact solution of the GP equation. Dashed lines: 
Gaussian approximation. 



FIG. 6: Stationary solutions of the GP equation versus the 
particle number J\T for a non-isotropic potential case with 
io r — Co/5 and lo z = Co (pancake). Top: value of the en- 
ergy functional. Bottom: square of the bifurcating eigenvalue 
(Aj_). Solid lines: exact solution of the GP equation. Dashed 
lines: Gaussian approximation. 



merical data. Comparing both results, we find that the 
Gaussian approximation captures the bifurcation quali- 
tatively, but with quantitative errors of 17% for A/~ r 1141 . 
24% for £a and 14% for A^ in the isotropic case (la ]. 

The phase portrait of the normal form is shown on Fig- 
ure|H When 5 = (1 -N/N c ) > 0, equation J2HJ) admits 
two fixed points Q± — T\/S/P, as shown in Fig. UJV. 
Thus, a hyperbolic stationary state and an elliptic sta- 
tionary state coexist. The phase space is separated into 
two regions by a separatrix which is a homoclinic orbit 
linking the hyperbolic stationary state to itself. Trajec- 
tories inside the orbit remain bounded near the elliptic 
fixed point. If the condensate is taken beyond the sep- 
aratrix by a perturbation (e.g. thermal excitations or 
quantum tunneling, see below section lTvjl . it will fall into 
unbounded (hyperbolic) trajectories and collapse. As j\f 
is increased, the hyperbolic and elliptic stationary states 
approach one another, (Fig. 0)3) and the homoclinic or- 
bit inside which orbits are bounded is reduced. The two 
stationary states join at Af — Af c (Fig.Qp), at which the 
HSN occurs. No stationary state exists for j\f > J\f c . 



The trial function is a Gaussian solution to the linear 
(a = 0) Schrddinger equation in which we incorporate 
eight variational parameters in order to take into account 
the anisotropy of the system. The form of the ansatz is 
given by: 

$(x,y,z,t) = (A r (t) + iAi(t)) 

where the real parameters {A r ,Ai}, {(j)x,4'Y,4>z} and 
{X, Y, Z} are related to the amplitude, the phase and the 
width of the Gaussian profile respectively. The Eulcr- 
Lagrange equations associated with the trial function 
(|35|l and the action defined in equation Q can be re- 
duced to the following system of second-order differential 
equations: 



d 2 X ~ v 1 

X -4- 

C. Non-isotropic Condensates dt 2 ~ x X 2 YZ X 3 

d 2 Y 2 _ v 1 



-ojf,Y - 



We now briefly present the main expressions obtained dt 2 v XY 2 Z Y 3 

from a Gaussian variational analysis of the G-P equation d 2 Z v \ 

with a cylindrical potential trap. Some of these results = ~ UJ zZ ~ z 2 XY 'Z 3 

have been prev iously obtained by other authors 0, ITsl 

where 



I30l I31L |32j, |33j • We therefore restrict our discussion to 
the equations that will be used in our analysis of the 1 2 m (uj x ujyUj z ) 1 ^ 3 

condensate lifetimes. v = V ~ t. \a\Af. (37) 
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The evolution of the condensate is better understood by 
drawing an analogy between its width and the motion 
of a particle with coordinates (X, Y, Z) moving in the 
potential 



TABLE I: Critical number of particles obtained for the 
isotropic, cigar and pancake geometries by using the exact 



solution of the GP equation (Af c 
imation (A/" C G ) 



and the Gaussian approx- 



U(X,Y,Z) 



1 
2 

2 { X 2 



ulX 2 
1 



i 

y2 



1 

z 2 



Indeed, defining P x = dX/dt, P y = dY/dt, P z = dZ/dt 
and the Hamiltonian 



V 


U r 


LJ Z 


Mi 


K G 


XYZ 


ill 


LJ 


1258.5 


1467.7 


(38) 


LJ 


co/5 


1460.3 


1646.6 




L>/5 


LJ 


1885.6 


2080.5 



H(P x ,P y ,P z ,X,Y,Z) 



1 



{Pi 



P^ + P 2 Z )+U(X,Y, 



we find that equations Ij36|) transform into Hamilton's 
equations of motion. 

If we consider now a potential trap JHJ with cylindrical 
symmetry (u> r = lo x = u> y ) equations i|36|l can be simpli- 
fied by using X(t) =Y(t). We thus find that ® yields 
two fixed points (X + ,Z + ) and (X_,i>_) which describe 
the stationary solutions for ^(x, y, z, t). These obey 

2 V 5 , : 



TABLE II: Rescaling factors required for having Afc = 
r?\^ in the cigar and pancake cases. Columns £ a and A A show 
critical amplitudes at the bifurcation for the rescaled £ and 
A 2 curves respectively. 







ixl 



n 2 4 M 

= U), — 



7Z1 



7Xi 



7Zi 



1 

7 \Xj 



1 



1 
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(39) 



where the chemical potential /i is related to the total 
number of particles through 



(40) 



The fixed points correspond to a metastable center 
(X+, Z + ) and to an unstable saddle point (X-, Z-), re- 
spectively. They are analogous to the Q+ an d Q- points 
appearing in the phase portraits on Figure 0] The so- 
lutions to (|39(l can be computed numerically, together 
with the linearized variational equations evaluated at ev- 
ery stationary point. 

Figures El and show £ and A 2 for the cigar and pan- 
cake cases, respectively. The solid lines present the val- 
ues obtained by discretizing and solving numerically the 
original differential equations (|12|l and (|22a(l . using the 
methods described in the Appendix. The dashed lines 
were computed using the Gaussian approximation de- 
scribed above. Both the isotropic and non-isotropic cases 
display saddle-node bifurcations. This is to be expected, 
since the saddle-node bifurcation is the gen eric way in 
which stable and unstable branches meet |23j. 

It is apparent from figures |3 03 and that the ex- 
act critical number of particles Aff is smaller than the 
Gaussian value A/" C G for all three geometries 0, 0,0,0] ■ 
Table|l]compares the different critical Af values obtained. 

In order to compare properly the HSN bifurcations ob- 
tained for the three aspect ratios studied, we can rescale 
the intensity of the potential to obtain the same A/" C E for 



LO T LO z 


rescaling factor £"a A a 


LJ LJ 


1340 14.68 


CcigW Ccigtll/5 


Ccig = 1.3463 1000 4.00 


Cpanti/5 Cp an [i 


Cpan = 2.2447 550 1.05 



all cases. In general, any confining harmonic potential 
with frequencies uj r and uj z that produces a critical num- 
ber of particles Af c can be rescaled by a factor 



Afc 



(41) 



to obtain the a new potential with frequencies lj* = c w T 
and uj* = coj z , which will have the critical number of 
particles Af*. The remaining physical quantities for the 
new potential are obtained through the following trans- 
formations: 



TV* 

£* 
X* 



7J74 
Af_ 

£_ 

X. 



(42) 
(43) 

(44) 
(45) 



We choose arbitrarily to rescale the potential intensity 
so that all Af® are equal to that for the isotropic case 
Nf a . Tablelrflshows the value of the rescaling factors c c - lg 
and Cp an for the cigar and pancake cases respectively, as 
obtained from equation glj using Af* = Af™ - The last 
two columns of this table show the critical amplitudes 
obtained for the rescaled £ and A 2 curves. These were 
obtained by fitting the HSN asymptotic forms given in 
relations l|3U|l and (|31|l to the rescaled data. 



IV. LIFETIME OF CONDENSATES 

In this section we first find expressions for the TIC, 
MQT and ICO decay rates. Using the numerical data 
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q 
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FIG. 7: Bounce trajectory (dashed) over the Euclidean po- 
tential <fr(q). Points g/, q m and (ft, indicate the fixed point, 
the minimum of <£>(?) and the bounce point, respectively. 



presented in the previous section, we then compute these 
decay rates for the u> r = uj z (isotropic), u r /5 — u) z (cigar) 
and u> r = lv z /5 (pancake) cases. Finally, we compare the 
results obtained for these three potential geometries by 
studying the spontaneous isotropization of the conden- 
sates. 



A. Definition and Computations of Decay Rates 

The TIC (thermally induced collapse) decay rate 
is estimated using the formula j2^| 



El 



lA-i 



2n 



exp 



—full 



(£+-£-) 



(46) 



where Huj(£ + — £—) is the (dimensionalizcd) height of the 
nucleation energy barrier (with uj the reference frequency 
introduced in Section III A|) , T is the temperature of the 
condensate and fcs is the Boltzmann constant. Note that 
the prefactor characterizes the typical decay time which 
is controlled by the slowest part of the nucleation dynam- 
ics: the top-of-the-barrier saddle point eigenvalue A+ and 
not A_ as used in [9j. However, near the bifurcation both 
eigenvalues scale in the same way and the behavior of Tt 
can be obtained directly from the universal saddle-node 
scaling laws (|30fl and (|31[l . Thus the exponential factor 
and the prefactor vanish respectively as 8 3 / 2 and 5 1 / 4 . 

We estimate the MQT (macroscopic quantum tunnel- 
ing) decay rate using an instanton technique that takes 
into account the semi classical trajectory giving the dom- 
inant contribution to the quantum action path integral 
8, 9]. We approximate this so-called bounce trajectory 
by the solution of the equation of motion 



d 2 q{t) _ -d$(q) 



dt 2 



dq 



(47) 



starting and ending at the fixed point g/ of the phase 
space where £(qf) = £-■ The Euclidean potential <t(q) 



is defined so that — "£(<?) reconstructs the Hamiltonian 
dynamics in the region scanned by the bounce trajec- 
tory (see Figure 0. We represent it by a fourth-order 
polynomial of the form 



&(q) = a + a 2 q 2 + a^q 3 + a 4 q 
coefficients ao, «2, c*3 and «4 chosen such that 



(48) 



*(0) = -£+ 

= -£_ 

d 2 m = -a+cao 



(49a) 
(49b) 
(49c) 

d 2 ®(q f ) = -A_(A0 (49d) 

We thus obtain a semi-analytic polynomial expression for 
$(</) where the coefficients are determined through the 
numerical values presented in figures |21 and |SJ 

Once $(?) and the bounce point qb (defined through 
the relation <&(<#,) = $(qf)) are known, the MQT rate is 
estimated as 



Co 



47T 



exp 



V2 



*(?) - &(q f )dq 



IS 



(50) 



where vq is defined by the asymptotic form of the bounce 
trajectory q(t) as it approaches g/ |9j, given by q(r) ~ 

9/ + (V|A-|)exp[-|A_T|]. 

In the same way as was done for the TIC, universal 
scaling laws can be derived close to criticality from 12811 . 
(|30|1 and i|31fl . We find that the exponential factor in 15U|) 
follows the same scaling as \/\£+ — £- \dq. It therefore 

vanishes as V 5 3 / 2 8 1 / 2 — <5 5 / 4 . The asymptotic form of 
q{t) shows that dq follows the same law as vo/|A_|. Thus 
vq ~ 5 3 / 4 , and the prefactor vanishes as V 5 1 / i 5 3 l i = 
<5 7 / 8 . Note that these universal scaling laws agree with 
those already derived in the Gaussian case in [g. 

The TIC (@EJ| and MQT (JH3J decay rates obtained for 
the exact and Gaussian stationary states are shown in 
Figure |H1 To validate these results we checked that the 
Gaussian TIC decay rates computed in 0] are found 
when we (incorrectly at a finite distance from criticality) 
replace A+ by A_ in equation l|46(l (data not shown). We 
also checked that our Gaussian MQT decay rate agrees 
with the one previously computed in p|. 

The ICO (inelastic two and three body collision) 
atomic decay rates are evaluated using the formula 
dM/dt = f c {N) with 



f c {M)=K / |*| 4 d 3 x + L 



(51) 



where K = 3.8 x Kr 4 [s~ 
[ill ITU . In order to compare the particle decay rate lj!)T|) 
to the condensate collective decay rates obtained for the 
TIC and MQT, we compute the condensate ICO half-life 



using 



ti/ 2 (A0 = 



.v 



dn 



A/72 fc(n) 



(52) 
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FIG. 8: Condensate decay rates versus particle number for 
the isotropic potential uo r = uj z — cj (solid), and for the 
rescaled cigar potential ui r = c c i g <i, u) z = c c i g ii/5 (dashed) 
and pancake potential lu t — c pan d)/5, ui z = c pa nW (dotted). 
ICO: inelastic collisions. MQT: macroscopic quantum tun- 
neling. TIC: thermally induced collapse at temperatures 50 
nK (3), 100 nK (4), 200 nK (5), 300 nK (6), and 400 nK (7). 




FIG. 9: Enlargement of the crossover region between the 
quantum tunneling and the thermal decay rate. ICO: inelas- 
tic collisions. MQT: macroscopic quantum tunneling. TIC: 
thermally induced collapse at temperatures 1 nK (1), 2 nK 
(2) and 50 nK (3). 

and plot Ty 2 on figure 00 

Figures 00 and 03 compare the condensate decay rates 
for the isotropic and the cigar and pancake potentials, 
rescaled by c c i g and c pan as described in Section IIII CI 
We note that the three aspect ratios generate very sim- 
ilar results after rescaling. The relative magnitudes of 
the different decay rates - TIC, MQT, and ICO - are the 




-0.2 0.2 0.4 0.6 0.8 1 



FIG. 10: Ellipticity ratio I as a function of fj,. Solid curves 
show numerical results, dashed curves the Gaussian approxi- 
mation. For the cigar, I = l r jt z (upper curves) and for the 
pancake, £ = £ z /l r (lower curves). Dots show results obtained 
by using successively fewer Fourier modes in the numerical re- 
sults; dots nearer to (further from) each curve correspond to 
retaining 7/8 (6/8) of the Fourier modes. £ changes by less 
than 1% for [i > —0.8 (fj, > —0.2) for the cigar (pancake) 
case and by less than 0.1% at the saddle-node bifurcation at 
fj, = 0.38 (p — 0.31) for the cigar (pancake) case. 



same for the three cases. At T < l[nK] the MQT effect 
becomes important compared to the ICO decay in a re- 
gion very close to (5 < 8 x 10 -3 ). This was shown 
in 0] using Gaussian computations but evaluating them 
with the exact maximal number of condensed particles 
Aff . Figure 00 shows that even for temperatures as low 
as 2[nK], the TIC decay rate exceeds the MQT rate ex- 
cept in a region extremely close to J\f c {5 < 5 x 10~ 3 ), 
where the condensates will live less than 10 _1 [s]. Thus, 
in the experimental case of 7 Li atoms, the relevant ef- 
fects are ICO and TIC, with the crossover that is shown 
on Figure 



B. Spontaneous Isotropization of Condensates 

The decay rates of the isotropic and non-isotropic cases 
shown in figure [S] are quite similar, despite the fact that 
lo z and lo t differ by a factor of five. We have investi- 
gated this question by examining the wave functions ^ 
for the pancake and cigar cases. These wave functions are 
peaked at the origin, as shown on figure 00 Their char- 
acteristic length scales in the axial and radial directions, 
i z and £ r , can be measured by computing the ratios of 
the value of ^ to its curvature at the origin. 
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FIG. 11: Length scales l T , lz as a function of [i. Solid curves 
show numerical results, dashed curves the Gaussian approx- 
imation. Cigar case shown in lower curves (with l r < l z ), 
pancake in upper curves (with l z < l r ). The size of the con- 
densate decreases drastically as /j, decreases, i.e. away from 
the saddle-node bifurcation along the unstable branch. 



More specifically, we define 

=)2, T .x -i 



02 - 



dz 2 



(r = 0,* = 0) 



(53a) 
(53b) 



We then obtain the ellipticity of the wavefunction as the 
ratio I of these length scales: t — £ r /£z for the cigar 
and £ = £ z /£ r for the pancake. These ellipticity ratios 
are shown in figure 1101 as a function of fx. For large fj,, 
i.e. away from the saddle-node bifurcation along the sta- 
ble branch, £ decreases rapidly away from one, indicating 
that the wave function is highly non-isotropic. At the 
saddle-node bifurcation, £ = 0.89 at \i = 0.38 for the 
cigar and £ — 0.80 at fx = 0.31 for the pancake. As fi 
is decreased, i.e. as we leave the saddle-node bifurca- 
tion along the unstable branch, £ approaches one as the 
wave function becomes more nearly spherically symmet- 
ric. This trend is present both in the numerical solution 
to the G-P equation and in the Gaussian approximation, 
as can be seen on figure [TJ3 Since the decay rates result 
from the scaling behavior near the saddle-node bifurca- 
tion, where the condensate is fairly isotropic, it follows 
that the decay rates are similar for the cigar, pancake, 
and spherically symmetric geometries, as we have shown 
on Figures |H1 and E| 

The spontaneous isotropization of condensates when 
fj, is decreased can be understood by the following phe- 
nomenological reasoning. When — /z grows, the balance 
of terms in the right hand side of equation JBJ changes. 
For small — fi, it is dominated by the the isotropic V and 



the anisotropic V(x) terms. But for large — //, the wave- 
function \I/ is strongly peaked and the V 2 and nonlinear 
terms, both isotropic, become dominant. 

Figure ITUI also provides a test of our numerical spatial 
resolution. By computing the ellipticity £ for different 
Fourier truncation levels, we show that £ changes with 
the resolution for low li, especially for the pancake case, 
where we used fewer Fourier modes than in the other cal- 
culations. Note however that our decay rate calculations 
only use results near the saddle-node bifurcation, where £ 
varies by less than 0.1% when different truncation levels 
are used. This indicates that 'F was adequately resolved 
in the region of interest. 

As [i is decreased, the wave functions become more 
highly peaked for both our numerical results and for the 
Gaussian approximation (see figures and 1111) . This is 
the main reason for the declining accuracy. To continue 
the computations further, the size of the periodic box 
should be reduced along with \i. We believe that, with 
adequate resolution, all of the exact wave functions would 
become spherically symmetric as fi decreases, as do the 
Gaussian approximations. 



V. CONCLUSION 

We have demonstrated that it is possible to numer- 
ically compute the stationary states, the bifurcating 
eigenvalues and the lifetime of anisotropic attractive Bose 
Einstein condensates. 

The Gaussian mean field approximation was found to 
have significant quantitative errors for all the different 
confining potential geometries that were studied, when 
compared with numerical solutions to the G-P equation. 

Spontaneous isotropization of the metastable conden- 
sate was found to occur as the critical number of particles 
is approached, yielding a life-time that depends weakly 
on the anisotropy of the confining potential. 

Direct methods - Gaussian elimination and diagonal- 
ization - were used in treating the spherically symmetric 
case, of size 128, but are far too costly for the three- 
dimensional case, of size 10 6 . In fact, since we only cal- 
culated axisymmetric stationary states and eigenvectors 
with an additional midplane symmetry, an intermediate 
two-dimensional axisymmetric cylindrical representation 
could have been implemented, of size 5000, permitting 
the use of direct methods. Our purpose, however, has 
been to construct and explore numerical methods appro- 
priate for a general non-isotropic case. 

The methods used to compute stationary states and 
bifurcating eigenvalues for the non-isotropic cases where 
essentially analogous. Each consists of a powerful and 
rapid outer iteration: Newton's method for the station- 
ary states and the inverse Arnoldi method for the eigen- 
values. The large linear systems that need to be inverted 
within each method are solved by the same inner bicon- 
jugate gradient iteration - BiCGSTAB - and constitutes 
the main numerical difficulty. Its convergence is greatly 
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improved by an inverse Laplacian preconditioning which 
is empirically tuned by adjusting the pseudo-timestep a 
in Newton's method or the shift s in Arnoldi's method. 

Our results and implementation have demonstrated 
that all these numerical techniques can be successfully 
combined to calculate the stationary states and eigenvec- 
tors for the G-P equation in a confining potential with 
an arbitrary three-dimensional geometry. 
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APPENDIX: NUMERICAL METHODS 
1. Spatial Discretization 

The operators L and W defined in l|10f) and l|ll(l are 

spatially discretized using the pseudospectral method 
|26j . For the isotropic case, the spherically symmet- 
ric ^(r, t) is expanded as a series of even Chebyshev 
polynomials T2 n (r/R), on which the boundary condi- 
tion ^(i?, t) — is imposed. The domain is taken 
tobeO<r<J? = 4 and the resolution used is 
Nr = 128. For the non-isotropic cases, we use a three- 
dimensional periodic Cartesian domain and '5 is ex- 
panded as a three-dimensional trigonometric (Fourier) 
series. The cigar case is solved in a periodic domain 
of size (L x , L y , L z ) — (5.39,5.39,12.04) in units of Lo, 
using (N x , N y , N z ) = (96,96,96) gridpoints or trigono- 
metric modes (with a 2/3 dealiasing rule,) so the total 
number of gridpoints or trigonometric functions is as high 
as Ns£> = 10 6 . (The more poorly resolved pancake case 
was calculated using (L Xl L y ,L z ) = (12.04,12.04,5.39) 
and (N x ,N y ,N z ) = (48,48,96).) The harmonic poten- 
tial (JSJ) is approximated by a periodic potential by writ- 
ing x = arcsin(sin(x)) and Taylor-expanding the arcsin 
function. This leads to a Fourier series for the potential, 
which is truncated according to the resolution used. 

Pseudospectral methods require performing over 4", 
at every iteration, a Chebyshev transform in the 
isotropic case or a Fourier transform in the non-isotropic 
case. These operations consume a time proportional to 
NftlogNji or N3D log(iV3r>), respectively. Actions and 
inversions of the Laplacian L are carried out on the 
Chebyshev or Fourier representations of 'J/, while actions 
of the multiplicative operator W are carried out on its 
grid representations. The time required by these opera- 
tions scales approximately linearly in Nr or N^o- 



2. Stationary States 

As stated in section III Bl the stationary states of J6j) 
that correspond to minima of £ at a given value of J\f 
can be obtained by integrating to relaxation the diffusion 
equation 

— = L* + (A.l) 

where the initial data ^(t = 0) has a total number of 
particles J\f and the value of the Lagrange multiplier jj, is 
fixed during the relaxation by the condition dAF/dt = 0. 

To integrate (|A.1|I a mixed implicit-explicit first-order 
time-stepping scheme is used: 

*(i + cr) = {I - aL)-\l + aW)^(t) (A.2) 

where I is the identity operator. The Helmholtz operator 
(/ — cL) _1 is easily inverted in the Chebyshev or Fourier 
representation. The motivation for integrating L implic- 
itly is to avoid the extremely small timesteps that would 
otherwise be necessitated by the wide range of eigenval- 
ues of the Laplacian. 

This relaxation method is equivalent to that used in 
[Tij and can only reach the stable stationary solutions 
of IjA.ljl . In order to also capture unstable stationary 
solutions EFil w e implemented a Newton branch-following 
algorithm |l5l |28| . We search for fixed points of l|A.2|) , a 
condition strictly equivalent to the stationarity of ©: 

Q = BV(t) = #(i + <r) - 

= (I -aL)- 1 ^ + aW)V{t) -*(<) 

= [{I -aL)- 1 ^ + aW)- I] #(t) 

= (I-aL)- 1 [(I + crW) - {I-<jL)]V(t) 

- (J-^^HL + Wp^). (A.3) 

Solutions to i|A.3|) are found using Newton's method. We 
begin with an initial estimate ^, in our case the solution 
at a neighboring value of /x. Newton's method calls for 
approximating the nonlinear operator B whose roots are 
sought by its linearization By about ^. We seek a decre- 
ment ip such that ^ — if) solves this linearized equation: 

= B($ - tp) w - B w tp 
B^iJj = B(V). (A.4) 

is then replaced by \& — -0 and equation l|A.4|) solved 
again for a further decrement. The process is iterated 
until B(i£) or ip is sufficiently small. In our case, equation 
(|A.4|I takes the form 

(I-aLY 1 a{L + DW)i\} = {I-aL)- 1 a{L + W)^. (A.5) 

We will explain how we solve the large linear problem 
(|A.5(1 in section HA 

The role of is formally that of the timestep in (|A.2(1 . 
but in HA.5fl . its value can be taken to be arbitrarily large. 
For a — > 00, (|A.5|I becomes: 

Lr l {L + DW)ip = L~ l {L + W)^. (A.6) 
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For the spherically symmetric case, the linear system 
HA. 611 is of size Nr — 128 and can be solved by standard 
Gaussian elimination. The boundary condition ip(r = 
R) = is imposed by modifying the operator L^ 1 or 
(I — crL)^ 1 , as it is in the time-stepping algorithm (|A.2|) . 

To compute the full branch of solutions as a function of 
/i, we begin from a stable state of l|A.l|) at a small value 
of M obtained by time-integration. Each stationary state 
is computed in 3 — 5 Newton iterations. 

3. Conjugate Gradient Solution of Linear Systems 

For the periodic Cartesian case, the linear system 
(|A.6|I of size N^£> = 10 6 is too large to be stored or 
inverted directly: the operation count for Gaussian elim- 
ination would be of the order of N^ D . Instead, we use 
BiCGSTAB [29). a variant of the well-known conjugate 
gradient method, developed for linear systems which are 
not symmetric definite. Such methods are matrix-free, 
meaning that they require only the right-hand-side of 
(|A.6|I , and a subroutine which acts with the linear opera- 
tor of the left-hand-side. A solution to the linear system 
is constructed as a carefully chosen linear combination of 
powers of the linear operator acting on the right-hand- 
side. 

For a periodic Cartesian geometry, conjugate gradient 
methods are particularly economical, since operator ac- 
tions are all accomplished in a time proportional to N^u . 
However, conjugate gradient methods for nonsymmetric 
definite systems may converge slowly (requiring a large 
number of evaluations of the linear operator) or even 
not at all. This happens when the operator is poorly 
conditioned, i.e. roughly when it has a wide range of 
eigenvalues. One must then precondition the linear sys- 
tem, i.e. multiply both sides of the system by a matrix 
which improves its conditioning and accelerates conver- 
gence. Since for operators such as L + DW, the wide 
range of eigenvalues is due primarily to those of L, we ex- 
pect L _1 to be an effective preconditioner. From i|A.5|) . it 
can be seen that a allows us to interpolate between linear 
operators a(L + DW) and L~ 1 (L + DW). We vary a em- 
pirically to optimize the convergence of BiCGSTAB. A 
few hundred BiCGSTAB iterations are usually required 
to solve the linear system. 

A further advantage of iterative inversion methods is 
that they can produce a (non-unique) solution even when 
the linear operator is singular. This is the case for our 
operators, which have the neutral modes described in 
section III CI as well as other neutral modes related to 
symmetries and the Fourier representation. The precon- 
ditioner, however, is inverted exactly. If l|A.6|) is used, 
the constant Fourier mode is treated separately which 
allows us to construct an invertible version of the singu- 
lar operator L. 



4. Eigenvalue Problem 

We now describe our numerical method for calculat- 
ing the linear stability of the stationary states. For the 
spherically symmetric case, the eigenvalues of (|17fl are 
computed by constructing and diagonalizing the corre- 
sponding matrix for each converged stationary solution. 
The results reported were generated with a Mathematica 
code running on a workstation. With the values R = 4, 
Nn = 128, the first two eigenvalues of the harmonic os- 
cillator are obtained with a precision better than 0.05%. 

For the three-dimensional case, it is again not possible 
to construct and diagonalize the matrix of size N^d = 10 6 
directly: the operation count for diagonalization is also 
of the order N^ D . Instead, we calculate only eigenvalues 
of interest, using a variant of the iterative inverse power 
method. The inverse power method calculates the eigen- 
values of a matrix M closest to a value s by means of the 
sequence defined by: 

(M - = 1>j (A.7) 

The sequence {tpj} converges rapidly to the eigenvector 
whose eigenvalue is nearest to s, with the eigenvalue A of 
M estimated by (A — s) _1 (tpj + x,tpj)/(tpj,tpj). 

In order to calculate complex or multiple eigenvalues, 
and to obtain more precise eigenvalues and error esti- 
mates, we use the sequence generated by (|A.7|) to imple- 
ment the more general Arnoldi or Krylov method [34j. 
Instead of retaining only the last two members of the 
sequence, the last K members (typically 4 or 6) are or- 
thonormalized and then assembled into the K xK matrix 
Hjk = (tpj, (M — s/) -1 ^)- The eigenvalues of H pro- 
vide estimates of up to K of the eigenvalues (A — s) _1 of 
(M-sI)- 1 . 

In our implementation of the Arnoldi method for (|22a|l . 
we seek the eigenvalues A 2 of the matrix — (L+DW I )(L+ 
DW ). Rather than solving 1A.7|) . we solve the equiva- 
lent preconditioned problem |28| 

^[-(L + DW^iL + DW^-sI]^^ =L- 2 ^ (A.8) 

by using BiCGSTAB. From l|A.8|) we obtain a sequence 
of vectors containing an increasing proportion of the de- 
sired eigenvectors, but since our solution of i|A.8|) is not 
exact, we then construct H by multiplication rather than 
inversion via Hjk = (tpj, Mtpk). We can then estimate 
the eigenvalues A 2 by those of H . Although the formal 
role of s is that of a shift which focuses the inverse itera- 
tion on the eigenvalues being sought, here we also use it 
empirically to improve the convergence of BiCGSTAB. 

The inverse Arnoldi method requires between 3 and 
10 iterations to converge, each of which requires several 
hundred BiCGSTAB iterations in order to solve its asso- 
ciated linear system. 
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